//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

// Polywell simulation. Simulate the movement of a single electron
// inside the polywell magnetic containment system.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#include "ephi.hpp"

// dihedral angle of the dodecahedron, 2arctan(φ), where φ=(1+sqrt(5)/2)
// as per http://en.wikipedia.org/wiki/Dodecahedron
#define DODECAHEDRON_DIHEDRAL_ANGLE 2.0344439357957027354455779231

void make_polywell_dodecahedron (Statics& statics, prec_t R, prec_t wire_radius,
    prec_t magnet_spacing, prec_t I, prec_t q_per_m, const CoilFactory& cf)
{
  // create a dodecahedron
  prec_t h2 = R * tan(DODECAHEDRON_DIHEDRAL_ANGLE/2);          // half of height
  vect3d mid = vect3d(0, 0, 0);                   // middle point of the dodecahedron
  prec_t x = R + R * cos(PREC_PI - DODECAHEDRON_DIHEDRAL_ANGLE);
  prec_t y = R * sin(PREC_PI - DODECAHEDRON_DIHEDRAL_ANGLE);
  vect3d bottom = mid - vect3d(0, 0, h2);         // bottom face center
  vect3d top = mid + vect3d(0, 0, h2);            // top face center

  prec_t real_r = R - (2 * wire_radius + magnet_spacing) * sin(DODECAHEDRON_DIHEDRAL_ANGLE/2);

  Ephi::debug (Ephi::DEBUG_INFO, "True magnet radius due to dodecahedron truncation and spacing: %f\n", prec2double(real_r));
  cf.make (statics, bottom, mid - bottom, real_r, wire_radius, I, q_per_m);
  for ( size_t i = 0; i < 5; i++ )
    {
      prec_t angle1 = (prec_t)i * PREC_PI * 2.0 / 5.0;
      vect3d bfc = bottom + vect3d(x * cos(angle1), x * sin(angle1), y);
      prec_t angle2 = angle1 + 2.0 * PREC_PI / 10.0; // offset from the angle1
      vect3d tfc = top + vect3d(x * cos(angle2), x * sin(angle2), -y);
      cf.make (statics, bfc, mid - bfc, real_r, wire_radius, I, q_per_m);
      cf.make (statics, tfc, mid - tfc, real_r, wire_radius, I, q_per_m);
    }
  cf.make (statics, top, mid - top, real_r, wire_radius, I, q_per_m);

  Ephi::debug (Ephi::DEBUG_INFO, "Polywell dodecahedron created, h = %f\n", prec2double(h2 * 2)); 
}

#define CUBE_DIHEDRAL_ANGLE (PREC_PI/2)

void make_polywell_cube (Statics& statics, prec_t R, prec_t wire_radius,
    prec_t magnet_spacing, prec_t I, prec_t q_per_m, const CoilFactory& cf)
{
  prec_t cubeR = R + (2 * wire_radius + magnet_spacing) / PREC_SQRT2;
  vect3d mid(0, 0, 0);
  vect3d pos;

  if ( !cf.make_coilpair (statics, vect3d(0, 0, 0), vect3d(0, 0, 1), 2*cubeR, R, wire_radius, I, q_per_m) ||
      !cf.make_coilpair (statics, vect3d(0, 0, 0), vect3d(0, 1, 0), 2*cubeR, R, wire_radius, I, q_per_m) ||
      !cf.make_coilpair (statics, vect3d(0, 0, 0), vect3d(1, 0, 0), 2*cubeR, R, wire_radius, I, q_per_m) )
    {
      Ephi::debug (Ephi::DEBUG_CRITICAL, "CoilFactory failed to make the coil pairs!\n");
      exit(-1);
    }

  Ephi::debug (Ephi::DEBUG_INFO, "Polywell cube created!\n");
}

void make_solenoid (Statics& statics, size_t n, prec_t R, prec_t wire_radius, prec_t spacing, prec_t I,
    const CoilFactory& cf)
{
  cf.make (statics, vect3d(0, -1e-7, 0), vect3d(1.0, 0.0, 0.0), R, wire_radius, I, 0);
  for ( size_t i = 1; i < n; i++ )
    {
      prec_t d = (prec_t)i * spacing;
      cf.make (statics, vect3d(-d, -1e-7, 0), vect3d(1.0, 0.0, 0.0), R, wire_radius, I, 0);
      cf.make (statics, vect3d(d, -1e-7, 0), vect3d(1.0, 0.0, 0.0), R, wire_radius, I, 0);
    }
}

